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Abstract 

A multigrid acceleration technique developed for solving the three-dimensional Navier- 
Stokes equations for subsonic/transonic flows has been extended to supersonic/hypersonic flows. 
An explicit multistage Runge-Kutta type of time-stepping scheme is used as the basic algorithm 
in conjunction with the multigrid scheme. Solutions have been obtained for a blunt conical 
frustum at Mach 6 to demonstrate the applicability of the multigrid scheme to high-speed flows. 
Computations have also been performed for a generic High-Speed Civil Transport configuration 
designed to cruise at Mach 3. These solutions demonstrate both the efficiency and accuracy of the 
present scheme for computing high-speed viscous flows over configurations of practical interest. 
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NAS1-18605 while the author was in residence at the Institute for Computer Applications in Science and Enginccring(ICASE), 
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Introduction 


During the last decade or so, significant progress has been made in the field of compu- 
tational fluid dynamics (CFD) to have an impact on the design and analysis of aerodynamic 
configurations. Solutions of the Euler (inviscid) equations for essentially complete aircraft con- 
figurations [1-4] and the solutions of the Navier-Stokes equations for high Reynolds-number, 
viscous, transonic flows over aircratt components are now available in the open literature [5—9]. 
It is noteworthy that most of the efficient numerical schemes for solving aerodynamic flows 
rely on multigrid acceleration technique [2,8,9] to enhance the convergence rate. The multigrid- 
based schemes have the desirable property that the number of iterations required to achieve a 
steady-state solution is nearly independent of the mesh size for a given class of problems. Thus 
one can achieve essentially grid-converged, steady-state solutions, even for the numerically de- 
manding problem of high Reynolds-number, transonic, viscous flow over realistic aerodynamic 
configurations with a reasonable amount of computer resources [9]. 

Despite the progress achieved in solving transonic flows, the development of CFD methods 
for supersonic/hypersonic flows seems to be lagging behind at the present time. With the current 
interest in high-speed vehicles such as the High-Speed Civil Transport (HSCT) and the National 
Aero-space Plane (NASP), it is imperative that efficient computational algorithms be developed 
for high-speed flow regimes. In this paper, we will discuss the progress made in this direction 
using a multigrid-based central-difference scheme and present sample results for problems of 
practical interest. 

Governing Equations and Numerical Method 

The basic equations under consideration here are the unsteady Navier-Stokes equations. 
These are specialized to a body-fitted coordinate system (£, ?/, Q, where £, i ], and C represent the 
streamwise, normal, and spanwise coordinates, respectively. The ?/ coordinate lines are nearly 
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orthogonal to the wing surface. Since the dominant viscous effects for high-Reynolds-number 
turbulent flows arise from viscous diffusion normal to the body surface, a thin-layer assumption is 
employed here by retaining only the viscous diffusion terms in the tj — direction. These equations 
can be written in the conservation law form as: 
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where the dependent variable vector U is given by the relation 
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In the eq. (1), F,G,G V , and H are the flux vectors, and J is the Jacobian of the transformation. 
The complete forms of these quantities are readily available in [9], 

A pseudo time-stepping scheme based on a Runge-Kutta scheme [10,11] is used for integrat- 
ing the time-dependent equations to steady state and as a smoother in the multigrid scheme. For 
convenience, let us first write the discretized form of the governing equations in the following 
operator notation: 


cl 

dt 




+ Q(U) - D(U) - 0 


( 3 ) 


where Q contains all the convective and viscous fluxes and D represents the artificial dissipative 
fluxes. 

Since our primary interest here is to obtain solutions for viscous flows via the Navier-Stokes 
equations, both the diffusion and the convective terms are important in contrast with the Euler 
equations where convective terms are dominant. Therefore, it is preferable to employ a scheme 
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that has a larger stability bound along the real axis in addition to good stability properties 
along the imaginary axis. Based on the Fourier stability analysis of a one-dimensional model 
problem, the five-stage Runge-Kutta scheme proposed by Jameson [12], with 3 evaluations of 
the dissipative operator at the first, third, and fifth stages, appears to be very attractive and is 
employed in the present work. The convergence to steady state is enhanced via the use of local 
time-stepping and implicit residual smoothing techniques [10,11], with the coefficients of the 
residual smoothing computed in the manner described by Vatsa and Wedan [9]. 

Artificial Dissipation Model 


The basic artificial dissipation model used in this study is patterned after the work of Jameson, 
Schmidt and Turkel [10] and of Jameson and Baker [11] for 2-D and 3-D Euler equations 
respectively, and modified later by Vatsa and Wedan [9] for 3-D Navier-Stokes equations. In 
order to discuss the modifications required ior supersonic/hypersonic flows, let us first examine 
the dissipation terms in the i-direction: 
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In the above expression, X i+ y <k is the modified eigenvalue scaling factor [9] and the coefficients 
e ( 2 ) and e ^ are related to the pressure gradient as follows: 


•+ jjit 
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(5) 


where the coefficients k.< 2 ) and k< 4 > are set equal to 1/2 and 1/64, respectively. The term «/ 
depends on the pressure gradient and is modified to give a TVD variation of the shock switch 
[13] in the following manner: 
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Note that by setting u> = 1, we can recover the shock switch that has been used in earlier studies 
for computing transonic flows [9-12], For supersonic and hypersonic flows, where shocks are 
much stronger, we use u = 1/2. The expressions for the dissipation terms in the j and k 
directions are derived in a similar manner. 

Evaluation of Time-Step 


It is very important to estimate the allowable time-step as accurately as possible in order 
to construct a robust time-stepping scheme. Failure to do so generally creates difficulties when 
the scheme is applied to different flow problems with widely varying test conditions and grid- 
densities. An attempt is made here to derive the expressions for allowable time-step for the 
present scheme from stability considerations. For convenience, let us start with the Navier- 
Stokes equations written in non-conservative form: 
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where U is the velocity vector and A, B, G zx are the coefficient matrices, their full form 

being available in [14 ]. 


Transforming these equations to the body-fitted curvilinear coordinate system £, r/ and (, 


and making the thin-layer assumption, one gets 


dU 

dt 


+ + B^y + + [Aij r + Br/y + Cy z }^~ 

( K oy 

+ [A( x + B( y + CQ~ 


‘Vi t y 

[D’ll + By 2 + Fyi + 6y v // r // v + G yz y y yz + G zx y z y x ] 


( 8 ) 


In general it is very difficult to derive an exact expression for the time-step At, since the 


coefficient matrices are non-symmetric. Abarbanel and Gottlieb [14] have recently delineated a 
procedure to symmetrize all of these coefficient matrices simultaneously. Taking advantage of 
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the work of [14], and the property that the norm of a symmetric matrix equals its spectral radius, 
one can find an upper bound on At in the following manner: 


— > — + + + (9) 

At A/c At,, A tc_ A tuff 

where, the first 3 terms on the right hand side arise due to the convective terms and the last term 
is due to diffusion terms. Setting the bounds of these components of At equal to their respective 
spectral radii, we arrive at the following expressions for the convective terms: 
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where c is the speed of sound. The diffusion limit on the time step is obtained in a similar 
manner and is given by the relation: 
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For viscous flow problems, the most restrictive time step is in the boundary layer region near 
a solid surface, where At^ and At^f are the most critical terms in determining the actual value 
of At. For adiabatic, high Mach number flow near a stagnation point, one can show that [13] 


A diff 


oc A 


(14) 


Thus the diffusion limit on time step becomes increasingly more important as the flow Mach 
number is increased. As expected, the diffusion limit reduces the allowable time step in the 
near wall region. However, by incorporating the diffusion limit in the time step, we were able 
to perform the computations over a large range of Mach number without changing the CFL 
number, which is very desirable from a practical point of view since it helps to make the scheme 
more robust. 
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Multigrid Acceleration Technique 


The convergence acceleration due to the use of multigrid techniques has been demonstrated 
for both inviscid and viscous flows [8,9, 12 j in the transonic flow regimes. In the current 
application, the Full Approximation Storage (FAS) scheme of Brandt [15] is used in conjunction 
with the multigrid strategy devised by Jameson [12] for the solution of the Euler equations. 
The extension of the scheme of [12] to the three-dimensional thin-layer Navier-Stokes equations 
for transonic flows described by Vatsa and Wedan [9] is used as a starting point for this work. 
A 5-stage Runge-Kutta scheme with coefficients selected to provide optimum damping of the 
high frequency errors is employed. The restriction operator used to transfer the solution to a 
coarser grid is a volume-weighted average of the eight surrounding cell-centered values. The 
forcing function for a cell on the coarse grid is obtained by simply summing the residuals of 
its constituent fine-grid cells. The corrections are transferred from the coarse grid back to the 
fine grid (or prolonged) by simple trilinear interpolation in computational space. On highly 
stretched or nonuniform grids, this prolongation operator can introduce high-frequency errors 
back to the fine grid, causing degradation of the convergence rate. To prevent this, the coarse- 
grid corrections were processed through an implicit residual smoothing operator before adding 
to the fine-grid corrections. Whereas the smoothing of the coarse-grid corrections was certainly 
helpful for transonic flow calculations, it was found to be essential for obtaining converged 
solutions for higher speed flows. 

The solutions presented in this paper were obtained using a W-cycle, in which governing 
equations are solved only in the restriction step. The W-cycle resulted in approximately a 
25-percent improvement in computational time compared to a standard V-cycle for achieving 
comparable convergence levels of the residuals. In addition, global properties such as lift and 
drag, develop more rapidly with the W-cycle, since more time is spent on the coarser grids. It 


6 


was also found helpful to run more cycles on coarser grid levels for supersonic and hypersonic 
flows in order to better establish and precondition the flowfield before starting computations on 
the finest mesh in the Full Multigrid (FMG) cycle. 

The variable-coefficient residual smoothings were applied on all grid levels of the multigrid 
cycle. On the finest grid, the blend of second- and fourth-difference artificial dissipation discussed 
previously was employed. For the coarser grids, a fixed coefficient second-difference dissipation 
model has been found adequate for transonic flow computations [9,12]. However, the coarse-grid 
dissipation model had to be modified via a pressure gradient based TVD switch (see eqs. (5-6)) to 
improve the convergence rate of the present scheme in supersonic and hypersonic flow regimes. 

Results and Discussion 

Two test cases covering supersonic to low hypersonic speed regimes are chosen for testing 
the multigrid Navier-Stokes code described in the preceding paragraphs. The accuracy of the 
computed solutions is assessed via comparisons with available experimental data. In the present 
investigation, C-0 type grids are employed. The computational grids are generated so as to 
cluster points in the appropriate regions to resolve sharp gradients present therein. In addition, 
significant grid clustering is used in the thin region adjacent to the solid surface in order to 
resolve the thin shear layers present in high Reynolds-number turbulent flows. 

a. Conical Frustum Entry Configuration 

The first test case used for this study is a simple aerodynamic shape designed for entry at 
Mach 6. The configuration consists of a modified conical-frustum [16] and is shown in Fig. 
(1), along with the downstream and symmetry planes. This configuration evolved through a 
conceptual study for the design of a vehicle to accommodate an 8— person crew, which could 
sustain a supersonic/hypersonic reentry and be capable of landing as a paraglider [16]. A wind 
tunnel model was tested at Moo — 6 0 and Reynolds number of 0.8xl0 6 based on the model 
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length. The Navier-Stokes calculations were performed on a grid consisting of 161x65x29 mesh 
points. The lift- and drag-coefficient data from the wind-tunnel tests are available for up to 
an angle of attack (a) of 12°. The Navier-Stokes solutions spanning this entire angle-of-attack 
range have been obtained to assess the performance of the current scheme over such a large 
range of test conditions. 

The convergence history in terms of the residual error of the continuity equation and the lift 
coefficient, C\ for the a = 6° case are shown in Fig. (2) as a function of work-units, where 
a work unit represents the computational effort required for one fine-mesh iteration. A total of 
400 iterations (620 work-units) were performed on the fine grid which resulted in approximately 
seven orders of reduction in the residual. The lift and drag coefficients for this case converged 
to within 0.1% of their final values in less than 50 fine-grid iterations. The convergence histories 
for this series of test cases are similar to the a = 6 U case shown in Fig. (2) except for the case 
of a = 12°, for which the residual started oscillating after dropping approximately five orders, 
possibly due to slight unsteadiness. 

The computed lift and drag coefficients for o = 0-12° shown in Fig.(3), are found to be in 
excellent agreement with the experimental data over the complete range of cv. Thus, not only does 
the present scheme have good convergence properties for these test conditions, but in addition it 
produces accurate solutions that are in good quantitative agreement with the experimental data. 

A better understanding of the overall flow field is obtained by examining Fig. (4), where 
the pressure contours for the a = 6° case are shown on the symmetry and downstream planes, 
in addition to the body surface. One can clearly observe the nearly conical growth of the shock 
surface in the streamwise direction. In the crossflow direction, the distance between the shock 
and the body surface goes through a minimum near the tip and then increases towards both the 
windward and the leeward planes of symmetry. 
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b. High-Speed Civil Transport 


The next test case considered here is that of the flow over a generic High-Speed Civil 
Transport (HSCT) configuration designed to cruise at a Mach number, M 0 0 = 3.0. The 
conceptual development and geometric details of this highly blended wing/body configuration 
are available in [17]. A wind-tunnel model representative of this vehicle was tested over a 
large angle-of-attack range for several Mach numbers, and the experimental data from this study 
have been documented by Covell et al. [18], In the present study, we will concentrate on 
the Moo = 3-0 case, which was the design cruise Mach number for this configuration. The 
corresponding Reynolds number based on model length was 6.3x10 . 

A grid consisting of 145x65x73 nodes was generated for the Navier-Stokes calculations. A 
partial view of the mesh on the symmetry plane and several streamwise cuts is displayed in 
Fig. (5). The outer boundaries of the grid were placed so as to contain the shock emanating 
from the leading-edge within the computational domain. Grid clustering in the tip, leading and 
trailing-edges are used to accurately resolve the flow in high gradient regions. 

The convergence histories for the lift and residual of the continuity equation for the a = 5 
are shown in Fig. (6), which shows a reduction of three orders in the residual after 300 fine 
grid iterations (470 work-units). The lift for this case converged after about 50 iterations. The 
convergence histories for the entire angle-of-attack range considered here ( a — 0 — b ) are very 
similar to this case. The integrated lift and drag coefficients are compared with the experimental 
data of [18] over the complete range of angle of attack in Fig. (7). It is clear from this 
figure that the computed results are in very good agreement with the experimental data. This is 
very encouraging especially for drag prediction, since viscous drag, which is difficult to predict 
accurately, constitutes a significant part of the total drag for this configuration at Mach 3. 

Next we examine the detailed surface pressure distributions for this vehicle. For this purpose, 


9 



we concentrate again on the c* = 5® case for which experimental pressure data are available at 
selected streamwise and spanwise stations. Since the computational grid does not follow the 
cuts along which experimental data was acquired, the computed pressure distributions had to 
be interpolated for a meaningful comparison. The interpolated solutions at two x-locations are 
compared with the experimental data in Fig. (8). For the station x=129, experimental pressures 
are available only on the upper surface, whereas both upper- and lower-surface pressures are 
available at the x=171 station. The agreement with the measured pressure data at both of these 
stations is quite good. It is noted that the computed solutions predict the correct variation 
in pressure even in the vortex-dominated flow near the wing-tip region. Similar comparisons 
have also been obtained at a = 1° and o = 3°; however these results are not shown here for 
conciseness. 

Whereas the global force coefficients and surface pressure comparisons are helpful in quan- 
titative validation of a prediction method, these flow properties are inadequate for understanding 
the true three-dimensional nature of the flow problem, such as the development of vortical-flow 
regions off the wing-tips. Following the lead of [19], an attempt is made here to visualize the 
vortical flow by plotting the density contours at fixed x-locations. This is done in Fig. (9) for 
x— 129 and x=171, the same stations for which the surface pressure distributions were exam- 
ined. The experimental laser sheet photographs, that are used routinely for visualizing shock 
and vortex formations are also shown in Fig. (9) at the corresponding stations. The core of 
the vortex forming between the wing-tip and the fuselage and its feeding sheet are the most 
dominant flow structure visible in these figures. After a careful examination of these results, 
the following observations can be made. First, the core of the vortex is predicted to grow in 
size and lift off farther away from the surface as one moves downstream. Secondly, the overall 
shape and size of the vortex predicted by the Navier-Stokes calculations are in reasonably good 
agreement with the laser-sheet data. 
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Concluding Remarks 


A multigrid acceleration technique developed originally for transonic flows has been extended 
to solve the three-dimensional Navier-Stokes equations for supersonic and hypersonic viscous 
flows. The convergence rate of the modified multigrid code for obtaining steady-state solutions 
of high Reynolds-number viscous flows in the supersonic/hypersonic flow regimes ( up to Mach 
6), is comparable to the convergence rate in transonic flows. Convergence histories and detailed 
comparisons with the experimental data are presented for two problems of practical interest. 
Based on these solutions, it is concluded that the resulting code is capable of predicting high- 
speed viscous flow problems in an efficient and reliable manner. 
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Comparison of force coefficient for conical frustum, 
.Uoc = 6.0 .Re I = O.SxlO 6 



17 





LOG(RESIDUAL) 



WORK 


Fig. 6 : Convergence history for IISCT on 145x65x73 grid. 
Moo = 3.0, « = 5 n ,Rci = 6.3x1 0 6 
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: Comparison of surface pressure distributions for HSCT, 
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